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In  this  report  we  disc 


:cts  of  the  method  of  domain  decomposition  for  the  solution 


of  elliptic  boundary  value  problems  in  two  dimensions.  The  basic  idea  behind  the  technique  is 
to  piece  together  local  solutions  of  the  elliptic  problem  to  form  a  global  solution.  There  are 
several  reasons  for  wanting  to  use  such  a  procedure.  The  first  is  that  a  given  domain  may  be 
irregular,  but  can  be  subdivided  into  regular  pieces  for  which  solutions  are  computationally 
efficient  to  obtain.  Another  reason  is  that  the  method  may  be  suitable  for  for  solving  elliptic 
problems  on  multiple  processor  machines.  If,  for  example,  one  breaks  up  the  domain  into 
many  small  pieces  and  allocates  a  processor  to  each  piece,  then  there  is  the  possibility  of 
decreasing  the  computational  time  by  constructing  the  local  solutions  in  parallel  and  "gluing* 
them  together.  \The  topic  of  domain  decomposition  has  become  increasingly  popular  recently 
and  many  papers  have  appeared  on  the  subject.  The  specific  domain  decomposition  method 


that  we  shall  be  concerned  with  is  that  discussed  and  analysed  by  Dryja  [8],  Bjorstad  and 
Widlund  [l],  Golub  and  Mayers  [9],  and  Bramble,  Pasciak  and  Hubbard  [2].  The  paper  by 
Bjorstad  and  Widlund  has  an  extensive  bibliography  as  well  as  a  short  summary  of  the  his¬ 
tory  of  the  method  The  purpose  of  this  report  is  to  examine  model  problems  with  the  inten¬ 
tion  of  providing  motivation  and  insight  concerning  results  previously  presented  as  well  as 


offer  some  new  observations  that  may  aid  in  the  implementation  of  the  method.  Ji - 
C  Earlier  papers  primarily  view  domain  deoompositionJ  from  a  discrete  viewpoint  (»  the 


method  is  applied  to  the  equations  resulting  from  a  discretisation  of  the  elliptic  problem. 
However,  underlying  the  method  for  the  discrete  equations  is  a  method  for  the  continuous 
equations.  We  shall  examine  the  method  for  the  continuous  equations.  We  are  interested  in 
this  viewpoint  because  it  provides  physical  interpretations  for  the  operators  that  arise  from 
the  discrete  application.  Also,  conclusions  derived  for  the  continuous  equation  should  be 
inherited  by  any  consistent  discretization.  For  problems  with  relatively  simple  geometry  (i.e. 
the  problems  we  are  considering  here)  explicit  representations  of  the  operators  occurring  in 
the  method  can  be  found.  This  aids  in  the  understanding  of  the  behavior  of  the  method, 
especially  with  respect  to  changes  in  the  geometry  of  the  domain. 


In  the  discrete  version  of  domain  decomposition,  one  is  faced  with  the  problem  of 
finding  the  solution  to  a  particular  set  of  linear  equations.  (The  equations  described  by  the 
capacitance  matrix.)  For  reasons  which  will  be  discussed,  it  is  advantageous  to  solve  the  sys¬ 
tem  using  an  iterative  technique.  The  iterative  technique  which  appears  to  be  the  most  popu¬ 
lar  is  the  method  of  pre-conditioned  conjugate  gradients  [4].  The  success  of  this  iterative 
scheme  relies  on  the  choice  of  a  pre-conditioner.  This  is  an  operator  whose  inverse  is  an 
approximate  inverse  for  the  linear  system.  It  is  our  aim  to  provide  some  insight  into  the 
nature  of  several  pre-conditioners  that  have  been  proposed,  and  the  effect  that  certain 
parameters  in  the  elliptic  problem  have  on  the  rate  of  convergence  of  the  resulting  iterative 
scheme.  Also,  perhaps  more  importantly,  we  shall  give  some  insight  into  the  physical 
significance  of  these  preconditioners  and  demonstrate  that,  although  they  may  be  defined  in 
different  forms,  they  are  all  essentially  the  same. 

The  continuous  approach  that  we  follow  uses  integral  equations.  We  also  work  with  the 
differential  form  of  the  equations.  Another  continuous  approach,  and  that  used  in  [l]  and  [2], 
is  the  variational  formulation  of  the  elliptic  problem.  We  use  the  integral  equation  approach 
because  in  the  case  of  simple  geometry  we  are  able  to  find  explicit  representations  of  the 
operators  involved.  We  also  found  it  more  useful  in  providing  motivation  for  the  method.  (In 
particular,  it  is  much  easier  to  interpret  various  preconditionere  using  such  an  analysis.)  The 
variational  formulation,  as  demonstrated  in  [l]  and  [2]  is  more  general,  being  able  to  deal  with 
variable  coefficient  elliptic  problems,  and  perhaps  more  useful  for  proving  certain  properties 
associated  with  method.  For  example,  Bjorstad  and  Widlund  [1]  prove  the  convergence  of  the 
conjugate  gradient  algorithm  for  the  capacitance  matrix  equations  when  their  suggested 
preconditioner  is  used.  It  is  possible  that  one  could  use  the  integral  equation  approach  given 
here  to  prove  results  concerning  domain  decomposition,  however,  it  seems  that  the  variational 
formulation  is  better  suited  to  such  a  task. 


In  the  first  section  we  give  a  brief  description  of  the  method  of  domain  decomposition 
from  the  discrete  viewpoint.  We  then  describe  the  method  from  a  continuous  viewpoint. 
Using  the  continuous  approach  we  examine  some  suggested  preconditioners.  The  continuous 
approach  is  then  used  to  examine  the  extension  of  the  domain  decomposition  ideas  for  a  vari¬ 
able  coefficient  elliptic  problem.  In  the  fifth  section  we  give  a  derivation  of  the  discrete  ver¬ 
sion  of  some  of  the  results  obtained  for  the  continuous  case.  This  discrete  version  is  based  on 
an  underlying  5-point  difference  scheme.  Pieces  of  this  derivation  first  appeared  in  [l]  and 
more  recently  a  complete  derivation  has  appeared  in  Chan  [6]  and  Chan  and  Resasco  [7],  We 
present  our  derivation  for  the  convenience  of  the  reader.  It  is  also  slightly  different  than  those 
previously  presented.  The  results  of  the  derivation  are  of  value  to  those  doing  domain  decom¬ 
position  on  rectangles  -  the  results  provide  the  information  necessary  for  an  implementation 
using  fast  Fourier  transforms. 
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this  report. 


1.  Description  of  the  Method 


In  this  Section  we  present  the  basic  ideas  of  domain  decomposition.  We  first  work  with 
the  discrete  version  of  the  method  and  then  the  continuous  version. 

The  problem  upon  which  we  shall  focus  primarily  is  the  solution  of  Poisson’s  equation 
on  a  domain  consisting  of  two  abutting  rectangular  regions.  Q  denotes  the  domain  and  0, 
and  O2  the  subdomains.  T  is  the  interface  between  the  two  subdomains.  (See  Figure  1.)  The 
equation  to  be  solved  is 

Au  =  / 

u  =  g  on  d  (1  (1.1) 

The  first  step  in  the  method  of  domain  decomposition  is  to  calculate  the  values  of  the  solution 
along  F.  The  two  subproblems 

Au,  «=  /  in  A,-  (1.2) 

u(-  =  u  on  T  and  u,-  =*  g  on  dfl,  /  T 

are  then  solved.  The  solution  on  fi  is  the  union  of  the  two  solutions  u,- . 

To  implement  this  idea  it  is  necessary  to  find  and  solve  an  equation  for  u  along  T.  We 
shall  call  this  equation  the  interface  equation.  It  may  seem  unlikely  that  such  an  equation  will 
exist.  What  we  shall  find,  at  least  for  the  model  problem  (and  it  can  be  expected  to  be  true 
for  other  problems),  is  that  there  is  an  equation  for  u  along  I\  The  form  of  this  equation 
depends  on  the  geometry  of  the  domain  while  the  influence  of  the  data,  /  and  g  ,  and  of  the 
values  of  u  in  each  of  the  subregions  will  enter  into  the  equation  as  forcing  functions. 

One  method  of  deriving  the  interface  equation  is  to  work  with  a  disc  re  tiled  version  of 
problem  (1.1).  This  approach  is  presented  by  Dryja  [8],  Golub  and  Mayers  [9],  and  Bjorstad 
and  Widlund  [1] .  Consider  a  standard  five-point  discretization  of  (1.1)  written  in  the  form 
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An  o  An 
0  Aj2  Agj 
An  Ajj 


Here  we  have  cancelled  one  factor  of  —  from  both  sides  of  the  equation.  The  values  of  u  are 

ft 

lumped  into  three  groups  Xj,  Xj  and  Xs  corresponding  to  points  in  Oj,  flj,  and  along  T 
respectively.  An  contains  the  coefficients  of  the  equations  for  the  solution  values  at  points  in 
Oi,  A-e  the  coefficients  for  the  values  at  points  in  ,  and  Ajj  contains  the  coefficients  for  the 
values  at  points  along  I\  The  coefficients  representing  the  coupling  between  the  different  sets 
of  points  is  contained  in  An,  Ajj  and  the  transpose  of  these  matrices.  One  derives  an  equation 
for  the  values  Xs  by  doing  block  Gaussian  elimination  on  the  matrix  system  (1.3).  (This  is 
just  Gaussian  elimination  on  (1.3)  thinking  of  the  sub- matrices  as  numbers.)  The  resulting 
equation  is 

[  C  j  j  =  y  £  Ass  -  AisA^Ais  -  A&AjaAgt  J  £  Xs  J  =*  Fs  -  A/sAjl'Fj  -  Aj»A^1F2 


This  system  of  equations,  represented  by  the  matrix  C,  is  the  Shur  compliment  of  Ass  with 
respect  to  An  and  A&.  C  is  often  called  the  capacitance  matrix.  It  can  be  shown  that  this 
system  is  negative  definite  if  the  original  system  is.  [5]  Equation  (1.4)  is  the  interface  equa¬ 
tion  that  we  are  seeking.  It  is  the  solution  of  this  system  which  is  the  principle  difficulty 
encountered  in  the  method  of  domain  decomposition.  Once  a  solution  is  found,  the  solution 
values  on  either  side  of  T  can  be  obtained  by  solving  the  equations  corresponding  to  each  of 
the  sub-domains  and  using  Xj  as  data  along  I\  Unfortunately,  the  matrix  C  is  difficult  to 
construct  explicitly  (in  general).  We  note  that  the  matrix  involves  Afi1  and  A® ,  which  are 
usually  only  known  implicitly  •  they  are  often  subroutine  calls.  Of  course  one  could  construct 
C  by  applying  C  to  a  set  of  basis  vectors  along  I\  Assuming  N  points  in  0  and  approxi¬ 
mately  N*  points  along  T,  this  would  involve  approximately  N*  applications  of  Aft1  and  A& . 


If  the  inverses  of  An  and  A&  were,  in  the  best  possible  case,  implemented  using  fast  solvers, 
then  the  resulting  work  would  be  the  order  of  Nl^aLogN.  This  may  be  acceptable  for  some 
problems,  but  the  attempt  is  made  to  solve  (1.4)  using  a  method  that  does  not  require  an 
explicit  representation  of  C.  In  light  of  this  requirement,  iterative  methods  are  a  natural 
choice  for  the  solution  of  (1.4). 

One  very  popular  iterative  technique  is  the  generalised,  or  preconditioned  conjugate 
gradient  algorithm  [4].  Given  the  problem  of  solving  Ax  =»  b,  one  can  think  of  the  general* 
ized  method  as  the  conjugate  gradient  scheme  applied  to  the  prepared  system 

M'1  Ax  =  M"lb 

in  the  inner  product  defined  by 

[*,  *]  =  (x,M-‘y) 

where  M  is  the  preconditioner  and  (-,  •)  is  the  standard  L2  inner  product  [2].  An  error  estimate 
[4]  for  the  resulting  iterative  procedure  is 


k  is  the  condition  number;  the  ratio  of  the  largest  eigenvalue  to  the  smallest  eigenvalue  of 
M*‘A 

Xa^pH^A) 

*  =  X^JM-'A) 

From  this  error  estimate  we  infer  that  a  good  preconditioner  is  one  which  is  a  multiple  of  an 
approximate  inverse  of  A. 

The  success  of  the  pre-conditioned  conjugate  gradient  method  is  dependent  upon  the 
choice  of  the  preconditioner.  For  the  model  problem  of  two  abutting  rectangular  regions, 
there  have  been  several  pre-conditioners  suggested  for  the  matrix  equation  (1.4).  We  shall  be 
interested  in  three  of  them.  The  first  two  assume  a  uniform  rectangular  grid  is  used  and  that 


the  nodes  on  the  interface  are  an  equal  distance  apart. 


Dryja’s  and  Golub  and  Mayers’  preconditioner  are  defined  in  terms  of  the  matrix  K, 
the  matrix  corresponding  to  the  standard  3-point  discretisation  of  the  one  dimensional  Lapla- 
cian  along  I\ 

-2  1 
1  -2  1 
1  • 

•  1 
1  -2 

The  suggested  preconditioner  is  to  use 

M  =  -{4Kf 

Here  (  )tt  represents  the  square  root  of  the  matrix  K.  The  second  preconditioner,  one  sug¬ 
gested  by  Golub  and  Mayers  in  [9],  is  a  modification  of  Dryja’s.  Specifically,  they  use 

M  =  -{4K  +  h  aK2)16 

The  third,  suggested  by  Bjorstad  and  Widlund,  is  to  use  as  M-1  the  operator  T  defined  impli¬ 
citly  as  follows: 

Let  Neumann  data  i?  be  specified  on  I\  Solve  in  Qj  the  problem 

86 

A  6  ”  0  with  on  T 

on 

6  — «  0  on  tin,  /  r 


and  setting  6  —  Ys.  Here  A^  is  the  appropriate  boundary  discretization  for  a  Neumann 
problem  on  flj.  There  is  nothing  special  about  using  flt  to  define  the  preconditioner,  one  could 
define  it  using  0$  as  well. 

For  convenience,  we  shall  refer  to  Dryja’s  preconditioner  as  D,  that  due  to  Golub  and 
Mayers  as  GM,  and  the  two  suggested  by  Bjorstad  and  Widlund,  one  defined  using  the  Neu¬ 
mann  -  Dirichlet  mapping  on  Qi  and  the  other  on  Cl?,  as  BWl  and  BW2  respectively. 

There  are  several  questions  related  to  the  use  of  these  preconditioners.  The  first  ques¬ 
tion  which  comes  to  mind  is  ”Why  are  these  good  preconditioners?*  Other  questions  concern 
the  relative  merits  of  the  preconditioners.  "Are  there  substantial  differences  between  them?” 
”Do  they  work  well  for  some  configurations  of  two  rectangles,  but  not  others?”  "How  easy  are 
they  to  implement?”  Also,  ”How  well  do  these  preconditioners  work  as  the  mesh  correspond¬ 
ing  to  a  discretization  tends  to  zero?”  And  lastly,  ”For  Bjorstad  and  Widlund’s  precondi¬ 
tioner,  which  piece  of  the  domain  should  be  used?” 

Some  questions  have  been  answered.  For  example,  in  the  case  of  the  problem  con¬ 
sidered  here,  numerical  experiments  presented  in  [9]  and  [1]  indicate  that  each  of  the  above 
preconditioners  is  very  good.  Typically,  the  conjugate  gradient  iteration  converges  in  a  lew 
iterations.  There  is  no  clear  best  preconditioner.  However,  it  appears  that  preconditioners 
GM  and  BWl  or  BW2  are  better  than  D.  As  for  theoretical  results,  Dryja  [8]  and  Bjorstad 
and  Widlund  [l]  have  been  able  to  obtain  bounds  on  the  condition  number  of  M_I  A  for 
preconditioners  D  and  BWl.  If  the  domain  given  in  Figure  1  is  discretized  using  a  standard 
five  point  approximation  to  the  Laplacian,  then  the  bounds  derived  by  Bjorstad  and  Widlund 
are 

+'t~+T^  +  2^2  U-6, 

k(BW1)  <  1  +  s/2  +  -(-?—■ - ) 
v  ?r  n  +1 

Here  p,q,r,m,  and  n  are  the  number  of  interior  nodes  corresponding  to  sides  of  the  domain.  In 
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Figure  2  the  correspondence  is  made  explicit.  Dryja  also  presents  a  bound  for  /c(D)  which  is 
similar  to  that  given  above.  One  interesting  aspect  of  these  bounds  is  the  condition  number 

for  both  preconditioners  becomes  unbounded  if  the  ratio  —  1  does.  Bjorstad  and  Widlund 

n  +  1 

conjecture  that  more  refined  estimates  may  not  have  this  behavior. 

For  problems  on  more  general  regions,  Bjorstad  and  Widlund  [1]  prove  that  for  the 
continuous  analogs  of  preconditioner  D  and  BWl  the  condition  number  of  M-1  A  is  bounded. 
(The  continuous  analog  used  is  based  on  the  equivalent  variational  formulation  of  the  prob¬ 
lem.)  If  the  bounds  (1.6)  are  sharp,  then  this  result  is  consistent  as  long  as  the  refinement  of 

the  underlying  mesh  is  carried  Out  in  such  a  way  that  —  +  —  and  —  remain  constant. 

n  +  1  r  +  1 

(This  is  a  very  natural  way  to  refine.) 

As  for  the  simplicity  of  computation,  if  one  has  a  problem  of  two  abutting  rectangles  in 
which  the  interface  is  a  straight  line,  then,  for  a  five  point  discretization,  all  of  the  precondi¬ 
tioners  can  be  implemented  using  fast  sine  transforms  on  T.  In  fact,  every  step  after  the  first 
of  the  conjugate  gradient  iteration  can  be  carried  out  using  fast  sine  transforms  on  the  line. 
For  the  specific  details  see  [5],  [6j  or  Section  5.  When  the  boundary  is  curved  it  is  no  longer 
possible  to  use  the  fast  sine  transforms.  Bjorstad  and  Widlund ’s  preconditioner,  not  depend¬ 
ing  upon  the  boundary  being  a  straight  line,  is  then  perhaps  the  easiest  to  apply.  All  that  is 
needed  to  implement  their  preconditioner  is  a  Neumann  solver  for  one  of  the  domains. 

To  help  answer  the  remaining  questions,  i.e.  motivation  for  the  preconditioners  and 
related  issues,  we  shall  consider  the  domain  decomposition  method  from  a  continuous 
viewpoint.  The  basic  procedure  is  the  same  as  in  the  discrete  case,  except  that  instead  of 
matrices  we  shall  have  integral  equations. 

For  a  domain  as  in  Figure  1,  the  principle  idea  is  to  construct  a  solution  which  is  the 
composed  of  two  solution  ut  and  U2.  These  functions  are  solutions  of  the  two  subproblems 
given  by  equation  (12).  The  boundary  data  for  Uj  and  u2  along  T  is  the  value  of  the  solution 
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u  along  T.  Over  the  rest  of  the  boundary,  U]  and  u2  take  the  boundary  values  g.  We  shall 
derive  an  equation  for  u  along  T  by  considering  the  conditions  that  ut  and  u2  must  satisfy  in 
order  to  form  a  solution  on  the  whole  domain.  Clearly  we  must  have  u;  and  u2  equal  along  T. 
The  other  condition  is  the  continuity  of  their  normal  derivatives.  This  condition  is  a  neces¬ 
sary  and  sufficient  condition  that  the  two  solutions  us  and  u2  must  satisfy  in  order  to  form  a 
weak  solution  of  the  equations. 

By  construction  Uj  and  u2  are  continuous  along  T.  What  is  not  guaranteed  is  that  the 
normal  derivatives  along  the  interface  will  be  continuous.  The  normal  derivatives  implicitly 
depend  on  the  values  of  u  along  T,  and  the  requirement  of  their  continuity  will  translate  into 
our  interface  equation.  We  now  use  Green’s  functions  to  find  an  explicit  representation  for 
this  equation. 

For  some  technical  reasons,  which  will  become  clear  later,  it  is  useful  for  us  to  consider 
as  our  unknown  the  difference  of  u  along  T  from  a  linear  function  which  interpolates  the 
boundary  values  of  u  at  the  end  points  of  T.  Thus,  we  shall  derive  an  equation  for  u'  where 

U*  =U  -  U| 

dtij  9u2 

We  derive  our  equation  by  finding  expressions  for  — —  and  — - —  along  T  and  then  set- 

on  on 

ting  them  equal.  We  denote  points  in  R2  as  x  =  (i  j,  x2)  or  y  =  (y  ,,  y2).  We  shall  consider 

d  d 

the  normal  derivative,  — — ,  to  be  taken  using  the  outward  normal  with  respect  to  ftj.  — - 

on  on, 

applied  to  a  function  of  two  sets  of  variables,  say  F(  x  ,  y ),  refers  to  a  normal  derivative 

du  j 

taken  with  respect  to  the  x  variables.  To  find  an  expression  for  — —  we  utilize  Green’s  iden- 

on 


Let  G]  be  the  Green’s  function  for  ft],  that  is 


A  Gj(z  ,  y)  =  6(x  -  y) 


■  W  W.5TK iKUTC  If* 


•  ■>  ■  •  j-  .-  . 


’  *.*  O  *.*  t.*  • 


1  •*  *  *  '  O'  v*  v*  v* 


and  G,(x  ,  y )  =  0  for  all  y  on  8Ctl.  Here  6  is  Dirac’s  delta  function.  We  have  the  identity 


dG 

Ui(* )  =  /  G,(x  ,  y )/  t(y  )dy  +  f  u  ,(s  )  (1.7) 

o,  so,  an » 

QQ  QQ 

=  /  G,(x ,  y )/  i(y  )<*y  +  /  si(« )  -r-1  (*,«)<*«  +  /  u'  (« (x , « )rf« 


/  !  is  the  restriction  of  /  to  fi,,  and  we  have  broken  up  the  boundary  values  of  u,  on  dOj 
into  two  pieces,  the  first  piece,  g,,  defined  by 


j  g  on  /  T 
Si  ~  (u|  on  T 

and  the  second  piece  equal  to  u'  on  I\  By  taking  the  normal  derivative  of  (1.7)  and  evaluate 
ing  it  along  T  we  have 


dui  8 

dnt  X  dnx 


dGx 


/ Gi(* ,  y )/  i(y  )dy  +  /  gi(< )  (* .  *  )* 


80, 


+ia' 


(We  have  assumed  that  it  is  valid  to  differentiate  under  the  integral  sign.)  A  similar  expres- 


chi2 

sion  holds  for  — —  in  f 1?.  Setting  these  expressions  equal  and  collecting  the  integrals  that 
on. 


contain  u'  on  the  left  hand  side,  we  obtain  an  integral  equation  along  T  for  u' 

|(x,  «)u'  («)dt  =  (1.8) 


<92G,  +  &Ga 


dn,  5S|  dn,  dn^ 


d 

dn. 


5G, 


/ G,(x  ,  y  )/  ,( y  )dy  +  /  g,(e  )  (x  ,  »  )dt 


80, 


dn. 


dG, 


/G  a(x ,  y )/  3(y  )dy  +  /  g 3(t )  -r—  (x,»)d» 

Oj  80j  dn» 


The  right  hand  side  of  (1.8)  has  the  interpretation  of  being  the  jump  across  T  of  the 


normal  derivative  of  the  solutions  to  two  Dirichlet  problems.  For  the  term 

/ Gi(* .  y )/  i(y  )dy  +  /  gi(« )  (x , • )d* 

Oj  an,  °nt 

is  the  normal  derivative  of  <f>,  the  solution  to 


a 4>  —  J  \  <t>  =  gi  on  an. 


(1.9) 


and  a  similar  interpretation  holds  for  Oj.  The  operator  on  the  left  hand  side  of  (1.8)  is  also 
composed  of  two  operators,  one  associated  with  each  domain.  The  operator  associated  with 

n,, 

j2q 

{ (’^  (iio> 

takes  data  u'  on  T  to  the  normal  derivative  of  the  solution  of 

A  V’,-  =  0  in  n<  (1.11) 

Vi  —  u'  on  T  Vi  =  Oon  O,-  /  T 

Since  each  operator  takes  data  on  T  to  flux  values,  we  shall  refer  to  these  operators  as 
’flux’  operators.  Thus  equation  (1.8)  expresses  the  fact  that  the  jump  of  the  normal  deriva¬ 
tives  (or  flux)  of  Uj  and  u2  induced  by  boundary  values  u'  along  T  must  equal  the  jump  in 
normal  derivative  (or  flux)  induced  by  the  forcing  function  fp-  in  each  region  and  the  boundary 
values  ~g  x  and  ~g  2. 

In  the  case  of  the  five-point  difference  approximation  to  the  Laplacian,  there  is  a  direct 


correspondence  between  the  matrix  problem  (1.4)  and  equation  (1.8): 


£  A53  ”  A|jA|ilA|3  “  A^A®1  Agj  j  * 

—  £  Am  +  »D  -  A/jAu'Ais  j  +  [  Am  +  »  D  -  A^A®  A#  J 

l  i 

{a3r<‘ '*<’>*  + 


and  for  the  right  hand  side, 


*  Fj  -  AjiAfi'F'i  +  hFj'  AjsAjaF2 


/G,(*  ,  1 )/  ,(>  )4t  +  /  gi(«  )  (*,•)*•  +  -g—  fG^x  ,  y  )f  j(y  )dy  +  /  g£» )  ~ 

a,  «i,  ani  an*  fij  «Oj  an 


Here  we  have  decomposed  AM  as  An  —  A«  +  D  +  A«.  Aj3  contains  the  entries  of  AM 
that  represent  the  coupling  of  the  nodes  of  T  to  0|.  Similarly  Am  corresponds  to  the  coupling 
of  nodes  of  T  to  Oj.  D  contains  the  entries  which  correspond  to  the  coupling  of  nodes  of  T  to 
each  other.  For  example,  consider  one  node  on  I\  then  the  entries  of  the  five  point  difference 
stencil  associated  with  the  node  would  be  broken  up  as 


2.  Explicit  Representation  of  the  Equations 

In  the  derivation  of  (1.8)  the  fact  that  the  domain  was  the  union  of  two  rectangles  was 
not  used.  However,  it  is  for  such  simple  regions  that  the  integral  equation  (1.8)  can  be  put  in 
more  explicit  form.  Before  we  consider  the  case  of  two  abutting  rectangles,  it  is  instructive  to 
consider  the  equation  for  an  even  simpler  geometry,  namely  that  of  the  half-plane. 

The  problem  to  be  solved  is 

in  =  /  in  R2 

We  consider  R2  broken  up  into  two  pieces  -  the  upper  and  lower  half  planes.  We  are  primarily 
interested  in  finding  a  representation  for  the  operator  on  the  left  hand  side  of  (1.8).  The  right 
hand  side  can  be  obtained  by  evaluating  the  jump  in  normal  derivative  of  solutions  of  a  Diri- 
chlet  problem  in  each  region.  (See  the  discussion  concerning  equations  (1.8)  and  (1.9).)  The 
Green’s  function  for  the  upper  half  plane  is  given  by 

G(X’  y)  =  ■j^Lo8([(*J-*i)  +  (*2-*2)]*  )--^L°8(I(*i-J'i)  +  (*8+  *2)1*  ) 


The  Green’s  function  for  the  lower  half  plane  is  similar.  If  one  follows  the  steps  that  led  to 
(1.8)  with  these  Green’s  functions,  then  one  arrives  at  an  integral  equation  of  the  form 


+OO  +00  ,  ... 

+  <2‘> 


where  R(s)  represents  the  jump  in  normal  derivative  along  the  x  j-axis  due  to  the  half-space 
problems 


A rl>i  =  /,  in  0f 
rp  =  0  on  x  j  =  0 

for  «  —  1,2.  n,  and  02  are  the  upper  and  lower  half  plane  respectively. 


The  interface  equation  can  therefore  be  represented  as  a  singular  integral  operator  equa¬ 
tion.  We  note  that  the  operator  (2.1)  is  a  convolution  with  a  shift  invariant  kernel,  and  hence 


the  Fourier  transform  of  the  equation  is  represented  by  multiplication  of  the  transform  of  « ' 
with  the  transform  of  the  kernel.  The  Fourier  transform  of  the  kernel  is  determined  by  the 
integral  operator  (2.1)  acting  on  functions  of  the  form  %'  (»  )  =  «*”** .  Instead  of  computing 
the  resulting  integrals  directly,  we  evaluate  the  integrals  using  the  knowledge  that  the  action 
of  the  operator  represented  by  (2.1)  can  be  computed  by  solving  the  two  problems  (1.11)  and 
evaluating  the  jump  in  the  normal  derivatives  across  the  x  ,-axis.  For  data  of  the  form  e  2ra* 
these  two  problems  can  be  solved  explicitly  using  separation  of  variables.  If  we  take  the  nor¬ 
mal  derivatives  of  the  solutions  we  find 


ea*a« 

(«  -o2 


4t 


Air  |  k  |  e8rtb 


Denoting  the  Fourier  transform  by  '  we  conclude  that  the  transform  of  equation  (2.1)  is 

-4ir|  k  |  u'  (*)  =  £(*)  (2.2) 

The  transform  of  the  1-D  Laplacian  is  -ta8*3  so  the  operator  for  the  interface  equation 
has  the  interpretation  of  being  -2  (-A,_D)'4.  Thus  the  motivation  for  preconditioners  D  and 
GM  is  that  they  are  both  finite  difference  approximations  to  -2  (-A,_D)'ti  and  hence  are  good 
approximations  to  the  exact  inverse  for  the  flux  operator  for  the  half-space.  Assuming  that 
the  operator  for  the  half-space  problem  is  an  approximation  to  that  for  rectangles,  we  under¬ 
stand  why  the  preconditioners  work  well  for  the  case  of  two  abutting  rectangles. 

BWl  and  BW2  are  preconditioners  is  defined  by  the  operator  which  takes  Dirichlet  data 
on  x  ,-axis  and  maps  it  (via  the  solution  of  either  half  space  problem)  to  Neumann  data  on 
the  x  ,-axis.  This  mapping  is  exactly  one  of  the  terms  on  the  left  hand  side  (2.1).  Thus, 
the  transform  of  such  a  preconditioner  is  just 

-2jt  |  k  | 

i.e.  one  half  the  transform  occurring  in  (2.2).  Hence,  for  a  half-plane,  any  discrete  approxi¬ 
mation  to  this  mapping  would  be  a  good  preconditioner.  In  fact,  for  this  special  problem,  D 
and  GM  can  be  considered  specific  realisations  of  Bjorstad  and  Widlund’s  preconditioner. 


16 


In  the  case  of  two  abutting  rectangles,  things  are  a  little  more  complicated  but  not 
much.  As  mentioned  earlier,  the  operator  in  the  interface  equation  (1.8)  is  the  sum  of  two 
"flux”  operators,  one  for  each  domain.  Consider  the  region  Oj  in  Figure  1.  We  are 
interested  in  a  representation  for  (1.10),  the  operator  that  takes  data  on  T  into  the  normal 
derivative  of  the  solution  of  Laplaces  equation  with  such  data  as  a  boundary  condition.  As  in 
the  half-space  case,  one  can  find  an  explicit  integral  equation  representation  of  this  operator 
using  the  Green’s  function  for  a  rectangle.  However,  by  manipulating  this  representation  (in 
a  long  and  tedious  calculation)  one  finds  that  the  eigenfunctions  of  the  operator  are  given  by 

k  ITS 

sin( - )  k  —  1,  2,  •  •  •  .  Once  the  eigenfunctions  are  known,  it  is  a  simple  matter  to  solve 

«i 

Laplaces  equation  with  such  data  using  separation  of  variables.  If  one  takes  normal  deriva¬ 
tives  of  these  solutions  and  evaluates  them  along  T,  one  finds  that  the  eigenvalues  correspond- 

k  ITS 

ing  to  the  operator  acting  on  sin( - )  are 

«i 

k  it 

A  similar  result  holds  for  the  operator  corresponding  to  flj.  From  this  information  we  can 
express  the  interface  equation  as 

[  ST'1  A„  ST0,  +  Tfa  ST/,1  A,  ST,,  E„,  ]  u'  («  )  *  R(.  )  (2.3) 

Here  ST0,  and  ST,,  represent  the  sine  transform  of  functions  defined  on  intervals  of  lengths 
»]  and  /?j  respectively,  i.e.  for  STa,  we  extend  a  function  defined  on  T  to  an  odd  function  on  a 
segment  of  twice  its  length  and  then  use  the  sine  series  of  the  extended  function.  It  is  here 
we  utilize  the  fact  that  we  are  representing  our  unknown  function  as  being  the  difference 
between  a  linear  interpolation  and  the  solution  values  on  the  interface.  This  representation 
guarantees  that  u'  will  be  continuous  when  extended  to  be  odd  periodic.  The  operator  E0, 
represents  extension  (by  zero)  of  a  function  defined  on  T  to  a  function  defined  on  the  top  side 


1  +  e 


-3rtait 


■v; 


• 


.  -  ■- .  V  -- 


of  O2.  T fa  represent*  the  operator  that  truncates  a  function  defined  on  the  top  side  of  fij  to  a 
function  defined  on  T.  and  Kf  are  diagonal  with 


k  jt  I  1  +  e 


(■Mi  ,*  —  — 


it  jt  I  1  +  e 


k  =1,2, 


We  see  that  the  eigenvalues  of  the  flux  operator  corresponding  to  rectangles  are  the 
eigenvalues  for  the  half  plane  multiplied  by  ratios  of  exponentials.  In  these  exponential  fac- 

Of 2  ^2 

tors,  the  character  of  the  domain  enters  in  via  the  aspect  ratio’s  —  and  — .  If  the  aspect 

»i  Pi 

ratios  are  large,  then  the  eigenvalues  are  essentially  those  corresponding  to  the  half  plane 
problem.  Also,  we  infer  that  effect  of  the  geometry  of  the  domain  is  most  pronounced  upon 
the  eigenvalues  associated  with  small  wave  numbers. 
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3.  Model  Problem 

We  shall  now  concentrate  on  the  case  when  both  rectangles  completely  share  a  common 
side  (or,  =  -  see  Figure  3.)  We  shall  refer  to  this  problem  as  the  "model”  problem.  In  this 
situation  and  and  both  are  the  identity  operator.  Equation  (2.3)  can  be  diagonalized; 

the  eigenfunctions  are  sin(  )  and  the  eigenvalues  are  given  by  (2.4).  Furthermore,  in  this 

ai 

basis,  the  preconditioners  described  in  Section  1  have  a  continuous  analog  whose  transfrom  is 
diagonal.  The  eigenvalues  corresponding  to  each  of  the  preconditioners  are 

41  0>, GM)  (3.1) 


kx  I  1  +  e 


(BWl) 


Here  k  =  1,  2,  •  •  •  .  BW2  has  a  continuous  analog  like  that  of  BWl,  but  with  &  replacing 
a2.  The  eigenvalues  of  the  continuous  analog  of  Dryja’s  and  Golub  and  Mayers’  precondition¬ 
ers  coincide  since  they  are  both  finite  difference  approximations  of  the  same  operator. 

Using  (3.1),  (3.2)  and  (2.3)  we  find  that  the  eigenvalues  of  M-1A  for  the  continuous  ana¬ 
log  of  the  preconditioners  are 


(D,  GM) 


-Siloj 

-3**#g 

1  + 

«  | 

1  +  e 

®i 

-8»t  og  T 

-s  «■* 

1  - 

e 

1  -  e 

-8 rkfa 
1+  e  01 

-2t*  #5 

1  -  e  ** 


(BWl) 


Certain  properties  of  the  preconditioners  can  be  inferred  from  these  expressions.  For  a 

Q  o  02 

decomposition  in  which  the  ratios  —  and  —  are  sufficiently  large  the  eigenvalues  depend 
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very  weakly  on  the  characteristics  of  the  domain.  The  eigenvalues  are  very  close  to  that  for 
the  half-plane  case.  We  expect  that  all  three  preconditioners,  each  of  which  is  an  approximar 
tion  to  the  exact  inverse  will  do  very  well. 


The  eigenvalues  for  D  and  GM  are  bounded  from  below  by  2,  and  the  maximum  eigen- 

($2  Oj 

value  tends  to  infinity  as  the  ratios  —  or  —  tend  to  sero.  Thus  the  condition  number  also 

a,  at 

tends  to  infinity  as  one  or  the  other  aspect  ratio  becomes  small.  We  conclude  that  the  condi¬ 
tion  number  associated  with  preconditioners  D  and  GM  cannot  be  bounded  independent  of 
the  domain.  Similar  conclusions  hold  for  BW1  and  BW2.  It  is  worth  noting  that  the  condi¬ 
tion  numbers  of  BWl  and  BW2  have  singular  behavior  with  respect  to  only  one  of  the  aspect 
ratios.  For  example,  if  fc  is  fixed,  then  for  all  a2  such  that  a2  <  ^  we  have  ic(BWl)  <  2. 
However,  if  a3  is  fixed,  then  for  all  <  a2  there  is  no  bound  for  a(BWl).  This  implies  that  if 
we  are  solving  a  problem  in  which  one  rectangle  has  small  aspect  ratio,  say  Ox,  then  we 
snould  use  the  Neumann  problem  for  0!  as  the  preconditioner  -  BWl.  The  condition  number 
of  BW2  will  be  correspondingly  large. 


C>2  ft 2 

If  we  have  a  domain  for  which  ( — )  <  ( — )  then  a  comparison  of  D  (and  also  GM)  to 

ai 


BWl  can  easily  be  made.  In  this  case,  the  minimum  value  of  the  eigenvalues  for  D  is  2  and 
the  minimum  eigenvalue  of  BWl  is  bounded  below  by  1.  Also,  every  eigenvalue  of  D  exceeds 
that  of  BWl  by 


-2**02 
1  +  e  a‘ 


-2*  *o9 


1  -  e 


(This  is  just  the  ratio  of  (3.3)  to  (3.4).)  From  these  considerations,  we  find  the  following  rela¬ 


tion  for  the  condition  numbers 


-1 


>  I 

k(BW1)  “  2 


1  +  e 


-a»*  a. 


1  -  e 


«2  . 


If  —  is  sufficiently  small  then  the  conditioner  number  associated  with  the  continuous  ana- 


logs  of  D  will  always  be  larger  than  that  for  BWl.  If  the  error  bound  (1.5)  is  a  reliable  indica¬ 
tor  of  performance,  we  would  expect  BWl  to  perform  slightly  better.  This  is  intuitively 
correct  as  BWl  takes  into  account  features  of  the  domain,  whereas  and  D  and  GM  do  not. 


In  conclusion,  we  see  that  all  the  preconditioners  can  be  expected  to  perform  rather  well 
since  each  contains  some  approximation  to  the  operator  -{-A1_d))6  which  is  the  relevant  for¬ 
ward  operator  of  the  domain  decomposition  problem  for  the  half  space.  The  characteristics  of 
the  domain  enter  in  weakly  as  long  as  one  side  is  not  significantly  smaller  than  the  other.  If 
one  side  has  a  small  aspect  ratio,  the  condition  number  of  D  and  GM  will  become  large  as  will 
one  of  BWl  or  BW2.  For  the  model  problem  at  least,  it  is  wise  to  use  the  rectangle  with  the 
smaller  aspect  ratio  to  define  the  preconditioner  suggested  by  Bjorstad  and  Widlund. 

In  in  the  discrete  case  (using  the  five  point  Lapiacian),  one  can  obtain  a  representation 
similar  to  (2.3).  This  shall  be  done  in  Section  5.  In  the  case  in  which  the  two  rectangles 
share  a  common  side,  the  discrete  operators  can  be  diagonalized,  and  the  interface  equation 
can  be  solved  directly  -  no  iteration  is  required.  One  may  then  wonder  about  the  usefulness  of 
the  above  analysis.  We  shall  demonstrate  by  numerical  experiments  in  Section  6  that  conclu¬ 
sions  based  upon  a  domain  composed  of  two  rectangles  sharing  a  common  side  carry  over 
more  general  domains.  An  analysis  similar  to  that  given  above  for  the  continuous  operators 
(a  comparison  of  preconditioners  etc.)  for  the  discrete  case,  has  been  carried  out  by  Chan  in 


4.  A  Variable  Coefficient  Problem 

We  now  discuss  the  application  of  domain  decomposition  techniques  to  a  variable 
coefficient  problem.  We  consider  a  model  problem  which  consists  of  solving 

V  «  V  *  =  f  «  =  f  on  dfl  (4.1) 

on  a  domain  0  as  in  Figure  1.  The  coefficient  e  has  a  value  e  t  in  fli  and  a  value  C2  in  flj. 
We  shall  dispense  with  a  discussion  of  the  method  of  domain  decomposition  for  a  discretiza¬ 
tion  of  (4.1)  and  only  discuss  the  method  from  a  continuous  viewpoint.  Again,  the  goal  is  to 
construct  a  solution  from  the  solution  of  two  subproblems 

e  ,  Aa  i  —  f  i  a  *=  g ,  on  dflj 
e  2A«  2  =  /  2  a  =  g  s  on  dCl^ 

The  boundary  values  g  i  and  j2  are  not  determined  along  T  and  an  equation  must  be  found 
for  these  undetermined  values.  This  equation  will  be  derived  by  considering  conditions  that 
the  two  solution  a ,  and  a  2  must  satisfy  in  order  to  form  a  solution  on  all  of  0.  First,  the 


solution  should  be  continuous  across  T, 


a  i  =  a  2  on  T 


and  secondly 


da i  da  2  .  . 

should  hold  across  I\  As  in  the  case  with  constant  coefficients,  these  requirements  are  neces¬ 
sary  and  sufficient  for  a ,  and  a  2  to  piece  together  to  form  a  weak  solution.  This  first  equation 
is  automatically  satisfied  by  our  choice  that  a  j  and  a  2  take  the  same  data  along  I\  It  is  from 
the  second  requirement  that  the  equation  for  the  interface  values  will  come.  Since  the 
coefficient  for  each  of  the  subproblems  is  constant,  we  can  use  Green’s  functions  for  each 
region  to  derive  expressions  for  the  normal  derivatives  of  a  i  and  a  2.  Using  such  expressions 
in  (4.2),  and  rearranging  the  resulting  integrals  in  a  way  similar  to  the  problem  with  constant 
coefficients,  we  arrive  at  an  integral  equation  of  the  form 
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Oj  »«!  on» 
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Oj  90j  * 


(4.3) 


Here  we  have  decomposed  the  solution  along  T  into  u'  and  a  linear  function  which  interpo¬ 
lates  the  boundary  values  at  the  endpoints  of  T.  ff  1,  ~02,  /  1  *nd  /  2  ate  defined  as  in  the  con¬ 
stant  coefficient  case.  Since  the  geometry  in  this  model  problem  is  so  simple,  we  can  use 
Fourier  analysis  to  express  this  equation  in  a  more  convenient  form.  The  eigenfunctions  of 

k  ir$ 

the  operator  on  the  left  hand  side  of  (4.3)  are  given  by  sin( - )  k  =  1,  2,  The 


operator  (4.3)  is  diagonalizable  by  the  sine  transform,  and,  denoting  the  sine  transform  of  a 
function  by  *  ,  it  takes  the  form 


' 

-2 v*  Oj 

-2*MS  ' 

e  Jr 

I  +  e 

e2k  ir 

1  +  «  ** 
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-2 wk  <*j 

-2»*tf2 

1  -  e  01 

1  -  e 
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*(*)  =  R(*) 


(4.4) 


In  a  way  analogous  to  the  constant  coefficient  case,  tl{k )  is  the  sine  transform  of  the  right 

d<t>i 

hand  side  of  (4.3)  and  can  be  evaluated  by  transforming  the  jump  across  T  of  c,-  for 

»  =  i)  2.  Here  fc  corresponds  to  solutions  of  Poisson’s  equations  in  Gi  and  flj.  Using  this 
formulation  we  shall  examine  a  preconditioner  which  is  an  extension  of  the  ideas  of  Bjorstad 
and  Widlund  and  is  suggested  by  Bramble,  Pasciak  and  Hubbard  in  [3].  As  in  the  constant 
coefficient  coefficient  case,  each  of  the  operators  occurring  on  the  left  hand  side  of  (4.3)  has 
the  interpretation  as  being  a  map  from  Dirichlet  data  on  T  to  Neumann  data  on  T .  The  idea 
is  to  use  one  of  the  operators  as  a  preconditioner.  The  inverse  of  the  preconditioner  (the 
operator  of  most  interest  to  us)  can  be  implemented  by  just  solving  a  Neumann  problem  on 
one  of  the  domains. 
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For  our  model  problem,  the  eigenvalues  of  such  a  preconditioner  are  given  by  one  or  the 
other  of  the  terms  on  the  right  hand  side  of  (4.4).  Let  Mi  and  M2  denote  preconditioners 
defined  in  terms  of  the  Neumann  problem  corresponding  to  0]  and  0$  respectively.  The 
eigenvalues  corresponding  to  the  matrix  M,rlA  are 
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As  in  the  constant  coefficient  case  we  see  that  the  geometry  of  the  domain  has  a  weak 
effect  as  long  as  the  aspect  ratios  are  large.  For  fixed  aspect  ratios  we  can  clearly  see  the  effect 
of  the  coefficients  on  the  condition  number.  In  particular,  if  e2«e1  then  the  condition 
number  associated  with  Mt  will  be  very  close  to  1.  (Similarly,  if  c2«c !  then  k(M2)  =  1.) 
This  suggests  that  one  should  use  the  domain  associated  with  the  larger  coefficient  to  define 
the  preconditioner.  If  one  uses  the  "wrong”  side  to  define  the  preconditioner,  then  we  expect 
a  larger  condition  number.  (This  will  of  course  depend  on  the  aspect  ratios  as  well.) 

The  ability  to  completely  diagonaliie  the  system  of  equations  (4.3)  using  Fourier 
analysis  carries  over  into  the  discrete  case  if  one  uses  the  standard  five  point  difference 
approximation  and  a  uniform  grid.  Thus,  the  discrete  equations  for  the  model  problem  con¬ 
sidered  here  can  be  solved  directly.  One  could  consider  the  solution  of  this  model  problem  as 
a  preconditioner  for  the  case  in  which  the  region  is  like  that  in  Figure  1.  (Just  extend  fii  into 
n2.)  We  denote  this  preconditioner  by  EV.  (For  extension  -  variable.)  Some  numerical  experi¬ 
ments  using  Ml,  M2,  and  EV  will  be  presented  in  the  next  Section. 

For  a  general  variable  coefficient  elliptic  problem,  Bramble,  Pasciak  and  Hubbard  sug¬ 
gest  using  a  preconditioned  conjugate  gradient  algorithm  in  which  the  preconditioner  is 
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defined  in  terms  of  a  variable  coefficient  “approximate”  equation  [3],  The  approximate  equa¬ 
tion  is  one  chosen  so  that  the  techniques  of  domain  decomposition  can  be  readily  applied. 
One  possible  approximate  problem  consists  of  decomposing  the  domain  into  subregions  and 
defining  the  coefficients  of  the  approximate  operator  to  be  constant  on  each  of  the  sub 
regions.  For  solving  this  approximate  equation  the  preconditioner  EV  may  be  useful. 


■  ’ ■ »  v  v  '.1  V-V,'.* V ■' 


25 


5.  Derivation  of  the  Discrete  Operators 

We  now  derive  the  discrete  eigenvalues  and  eigenvectors  corresponding  to  the  capaci¬ 
tance  matrix  C  in  (1-4).  For  alternative  derivations  see  [l],  [6]  or  [7].  As  discussed  in  Section 
2,  the  continuous  analog  of  this  matrix  is  an  operator  composed  of  two  parts.  Each  part  is 
associated  with  a  domain  and  represents  the  mapping,  defined  in  terms  of  a  solution  of 
Laplaces  equation,  of  Dirichlet  data  on  T  to  Neumann  data  (or  flux)  on  T.  It  is  the  same  in 
the  discrete  case,  and  therefore  it  suffices  to  find  the  eigenvalues  and  eigenvectors  of  the 
discrete  operator  for  one  rectangle.  Our  derivation  will  essentially  follow  the  procedure  for  the 
continuous  case  given  in  Section  2.  Consider  the  rectangle  in  Figure  4.  Given  data  on  T,  the 
flux  induced  at  a  point  (lht  0)  is  given  by 

2)~  [*„o  ~  2^1.0  +  -n,o  J  +  [Vi,i  -  ^,o  j  (5.1) 

applied  to,  <j>,  a  solution  of 

A*  <l>  —  0  on  H,  (5.2) 

<t>  —  9  on  T 

Here  we  use  A*  to  denote  the  standard  five  point  Laplacian.  h,  and  A,  denote  the  mesh 
widths  in  the  x  and  y  directions  respectively. 

We  recognize  the  second  term  of  (5.1)  as  being  the  finite  difference  analog  of  the  normal 
derivative.  The  first  term  is  due  to  the  flux  induced  by  values  along  T  and  does  not  appear  in 
the  continuous  case.  However,  the  inclusion  of  it  is  natural  if  one  considers  the  derivation  of 
the  five-point  difference  stencil  from  a  control  volume  point  of  view. 

ickih 

The  eigenvalues  for  the  operator  represented  by  this  mapping  are  sin( - — )  i.e.  the 

at 

continuous  case  eigenvalues  evaluated  at  the  nodes.  This  can  be  deduced  from  a  manipulation 
of  the  discrete  Green’s  function  for  the  rectangle,  or  from  matrix  considerations  ([1] ,[6),[7]). 

We  now  precede  to  find  the  eigenvalues  associated  with  such  eigenfunctions.  We  shall 
do  this  by  solving  (5.2)  with  eigenfunction  data,  and  then  apply  the  difference  operator  (5.1) 


to  the  result.  Our  solution  of  (5.2)  will  be  constructed  from  solutions  for  a  domain  which  has 
a2  =  oo.  These  solutions  will  be  combined  in  such  a  way  that  the  boundary  condition  along 
the  top  of  the  rectangle  will  be  satisfied. 

For  the  half-infinite  rectangle,  we  employ  a  discrete  version  of  separation  of  variables 
and  seek  solutions  of  the  form 
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(5.3) 


We  have  denoted  the  dependence  of  X  on  the  wavenumber  k  by  adding  the  subscript.  The 


results  for  sin  ( 


nklh 


)  will  be  obtained  by  combining  solutions  of  the  form  (5.5)  for  k  and 


-k  .  We  determine  X*  by  substituting  (5.3)  into  the  difference  equation  (5.2).  At  the  point 
( lhs  ,  mht )  we  find 
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If  we  cancel  the  terms  e  1  and  X4"  and  simplify,  we  arrive  at  an  equation  for  X4 

X*2  -  (2  +  A(~fnn2(  2q*  ))X*  +  1  •=  0 


(5.4) 


We  use  \k  corresponding  to  the  root  of  (5.4)  which  is  less  than  1  (to  get  a  solution  bounded 
at  oo).  Thus, 


The  solution  for  the  finite  rectangle  is  now  an  infinite  sum  of  solutions  of  the  form  (5.3)  with 
X*  defined  above, 
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Here  N  is  the  number  of  panels  along  the  side,  i.e.  — — .  Formula  (5.5)  is  derived  using  the 

method  of  images.  The  simplification  is  obtained  by  using  the  summation  formula  for  a 
geometric  series. 

The  difference  operator  (5.1)  is  now  applied  to  the  solution  (5.5)  along  T  and  this  yields 
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The  eigenvalue  At  is  also  an  eigenvalue  for  data  of  the  form  e 
be  an  eigenvalue  for  a  linear  combination  of  the  two,  i.e.  sin( 
Jfc  =  1,  2,  •  •  •  are  the  eigenvalues  that  we  seek. 


-s» Mkt 
2a 

nklh. 


and  therefore  will 
).  Thus  A*  with 


Using  these  eigenvalues  and  eigenvectors  it  is  a  simple  matter  to  compute  the  discrete 
operator  corresponding  to  (1.4).  If  the  domain  is  that  given  in  Figure  1,  then  we  use  the  for¬ 
mula  (2.3)  with  discrete  sine  transforms  replacing  the  continuous  sine  transforms.  We  have  a 
representation  of  the  matrix  equations  as 

|  A33  -  AijAf/Au  -  A23A22  A23  J  X8  =  (5.7) 

[  ST*-;  A,  STa,  +  T*  ST/;  A,  ST,,  Ea,  ]  X,  =  Fs  -  A/sAf/F,  -  A^A^F2 

The  operators  and  vectors  on  the  right  hand  side  of  (5.7)  are  those  defined  in  (1.4). 
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STa],  ST(>i  represent  the  discrete  sine  transform  of  functions  defined  on  intervals  of  lengths  qx 
and  0i  respectively.  The  operator  Eaf  represents  extension  (by  xero)  of  a  function  defined  on 
r  to  a  function  defined  on  the  top  side  of  Oj.  T*,  represents  the  operator  that  truncates  a 
function  defined  on  the  top  side  of  fij  to  a  function  defined  on  T.  A<,  is  diagonal  with  its 
(k  ,k  )th  entry  given  by  A*  in  (5.6).  A f  is  diagonal  and  has  entries  similar  to  those  of  A«  but 
with  0i  used  instead  of  at.  The  right  hand  side  of  the  equation  can  be  formed  by  evaluating 
the  jump  in  flux  induced  by  the  solution  of  two  Poisson  equations.  (See  the  discussion  con¬ 
cerning  equations  (1.8)  and  (1-9).)  In  the  iterative  solution  of  equations  (1.4)  it  is  only  neces¬ 
sary  to  compute  the  product  of  C  with  a  vector.  From  (5.7)  we  see  that  this  can  be  done 
completely  using  discrete  sine  transforms.  Also,  if  0i  =  a,  then  the  discrete  operator  can  be 
diagonalized  and  the  equation  (1.4)  can  be  solved  directly. 

For  the  discrete  half  space  case,  the  eigenvalues  are  the  same  as  that  for  the  rectangle 
except  that  N  =  oo.  If  A,  =  ht  =  h  ,  then  the  eigenvalues  A*  simplify  and  are  given  by 


=  -V4K+  A2 K2 

m 


Since  the  eigenvalues  for  the  full  half  space  operator  are  twice  those  of  (5.8),  we  see  that 
Golub  and  Mayers’  preconditioner  is  the  exact  inverse  for  the  discrete  half  space  problem.  For 
more  comparisons  of  the  discrete  preconditioners,  see  Chan  [6]. 
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6.  Numerical  Experiments 

We  now  discuss  the  resuits  of  some  numerical  experiments.  In  all  cases  the  discretira- 
tion  will  be  that  arising  from  a  five-point  difference  approximation  to  the  Laplacian.  The 
mesh  is  uniform  with  width  h  in  each  direction.  Since  we  are  using  a  regular  distribution  of 
points  and  the  interface  is  a  straight  line,  the  iterative  solution  of  equation  (1.4)  can  be  car¬ 
ried  out  using  fast  sine  transforms. 

In  the  first  experiment  we  worked  with  a  domain  as  in  in  Figure  1.  We  were  interested 
in  the  behavior  of  the  condition  number  of  M_1A  for  each  of  the  preconditioners  as  the  ratio 

ai 

—  was  varied,  i.e.  as  the  region  becomes  more  ”L”  shaped.  The  goal  is  to  see  if  the  conclu¬ 


sions  based  on  the  model  problem  were  valid  for  perturbed  domains.  We  examined  five 
different  pre-conditioners;  D,GM,  BWl,  BW2,  and  one  which  was  the  exact  solution  for  the 
domain  consisting  of  the  extension  of  the  narrower  rectangle  into  the  larger  rectangle.  We 
shall  refer  to  this  as  preconditioner  EC  (Extension-Constant).  Chan  in  [6]  has  also  recom¬ 
mended  this  as  a  preconditioner.  Using  the  results  of  the  previous  Section,  the  application  of 
all  of  the  preconditioners  can  be  completely  carried  out  using  fast  sine  transforms.  When 

ai 

—  =  1  then  we  are  solving  the  model  problem.  The  results  of  the  experiment  are  plotted  in 
Pi 

Figure  5.  From  these  results  we  see  that,  except  for  BW2,  the  condition  number  given  by  the 
model  problem  accurately  predicts  the  value  of  the  condition  number  over  a  wide  range  of  the 


ratio  — .  It  appears  that  the  condition  number  appears  bounded  with  respect  to  perturba- 
P\ 

tions  in  this  ratio.  Except  for  minor  deviations,  BWl  and  G  have  about  the  same  condition 
number  while  that  associated  with  D  is  slightly  larger.  BW2  appears  to  be  worse  than  all  the 
other  preconditioners.  The  difference  between  k(D)  and  k(G)  perhaps  lies  in  the  accuracy 
with  which  each  approximates  the  discrete  half-space  flux  operator.  In  this  experiment  we 
chose  the  component  rectangles  such  that  cr2  =  &.  Some  experiments  with  a2y^  &  were  per¬ 
formed,  and  it  was  found  that  the  results  did  not  differ  significantly.  One  of  these,  that  in 


which  —  =  .1  is  presented  in  Figure  6. 
n 

In  the  second  experiment  we  again  used  a  domain  like  that  in  Figure  1,  but  varied  the 
°s 

the  ratio  — .  As  discussed  in  the  Section  3,  we  expect  that  as  this  ratio  tends  to  tero  the 

<*i 

condition  number  associated  with  BW2,  D  and  GM  should  grow.  The  results  for  two  different 
width  fingers  is  presented  in  Figures  7  and  8.  For  a  finger  which  is  3/S  the  width  of  the  bot¬ 
tom  rectangle,  we  indeed  see  the  expected  growth  of  the  condition  number  as  the  aspect  ratio 
a2 

—  was  made  smaller.  *c(BWl)  grows,  but  does  not  appear  to  do  it  catastrophically  as  it  does 

«i 

with  D,  GM  and  BW2.  We  also  note  that  our  fifth  preconditioner  EC,  does  very  well  -  in 
fact,  the  condition  number  gets  smaller.  For  a  very  thin  finger,  one  which  is  1/5  the  width  of 
the  bottom  rectangle,  the  growth  of  the  condition  number  occurs,  but  it  is  not  as  marked  as 
the  previous  case.  For  very  thin  fingers,  the  preconditioners  appear  to  become  uniformly 
better. 

The  above  estimates  of  the  condition  number  were  primarily  made  to  assess  the  predic¬ 
tive  power  of  our  analysis  for  the  model  problem.  However,  the  condition  number  is  only 
used  in  an  error  bound,  and  it  may  not  necessarily  indicate  the  true  performance  of  the  itera¬ 
tive  method.  To  see  if  the  condition  number  is  indeed  a  reliable  indicator,  we  determined  the 
number  of  conjugate  gradient  iterations  that  were  required  to  reduce  the  increment  to  a  value 
less  than  C  A  2.  (By  increment  we  mean  the  amount  added  to  the  current  approximation  to 
obtain  the  next  one.)  This  stopping  criterion  was  chosen  since  we  were  interested  in  solving 
the  equations  to  an  accuracy  which  was  a  multiple  of  the  truncation  error  of  the  numerical 
difference  scheme.  We  choose  the  constant  s C  (somewhat  arbitrarily)  to  be  I0_1.  The 
domain  chosen  was  that  which  gave  rise  to  Figure  7.  These  results  are  presented  in  Table  1. 
As  the  condition  number  of  every  preconditioner  is  constant  for  a  wide  range  of  the  parame¬ 
ter  it  is  no  surprise  that  the  number  of  iterations  was  essentially  constant  as  well.  In  the 
cases  where  the  condition  number  was  large,  there  was  an  increase  in  the  number  of 


iterations.  (One  or  two  more  iterations.)  It  appears  that  EC  does  uniformly  well  over  the  com¬ 
plete  range  of  the  aspect  ratio.  The  results  for  GM  and  BWl  are  almost  identical  to  EC 
except  for  fairly  small  aspect  ratios.  D  and  BW2  are  a  bit  worse,  requiring  on  the  order  of  two 
more  iterations  than  EC.  Compared  to  the  results  using  no  preconditioning,  all  the  precondi¬ 
tioners  do  rather  well. 

In  the  last  experiment  we  investigated  the  effect  of  variable  coefficients.  The  domain 
was  a  region  as  in  Figure  1.  The  width  of  the  centered  finger  was  3/5  that  of  the  bottom. 
(The  results  for  non-centered  fingers  were  nearly  identical.)  Three  preconditioners  were  used. 
The  first  two  were  Ml  and  M2  of  Section  4;  preconditioners  defined  in  terms  of  Neumann 
problems  on  and  O?  respectively.  The  third  one,  EV  of  Section  4,  was  the  exact  solution 
operator  for  a  domain  which  consisted  of  extended  into  ft?  and  using  the  coefficients  of  flj 
for  the  operator  in  the  extension.  Due  to  the  simple  geometry,  each  of  these  preconditioners 
and  the  iterative  scheme  for  solving  the  capacitance  matrix  equations  was  implemented  using 
fast  sine  transforms.  Our  stopping  criterion  for  the  iterative  scheme  was  the  same  as  that  for 
the  constant  coefficient  case.  In  the  first  experiment  e ,  —=  100.0  and  ca  —  1.0.  In  Table  2  we 

a2 

present  the  number  of  iterations  verses  the  ratio  —  for  each  preconditioner  and  the  standard 

®i 

(unconditioned)  conjugate  gradients.  As  discussed  in  Section  4,  the  coefficients  for  this  prob¬ 
lem  can  be  expected  to  reduce  the  condition  number  of  Ml.  This  indeed  appears  to  be  the 
case,  as  the  number  of  iterations  it  takes,  3,  is  only  one  more  than  can  be  expected  with  the 
exact  inverse  as  a  preconditioner.  Our  suggested  preconditioner  EV  did  as  well,  while  using 
M2  was  particularly  bad.  The  results  for  M2  are  expected  since  the  aspect  ratio  as  well  as  the 
coefficients  have  an  unfavorable  effect  on  the  condition  number.  In  the  second  experiment, 
e  i  =  1.0  and  e2  — =  100.0.  The  results  are  presented  in  Table  3.  The  effect  of  the  coefficients 
can  be  seen  by  an  increase  in  the  number  of  iterations.  For  Ml,  depending  on  the  aspect 
ratio,  it  either  doubled  or  tripled  the  number  of  iterations.  This  clearly  suggests  that  if  one 
uses  a  preconditioner  defined  in  terms  of  a  Neumann  problem,  the  Neumann  problem  should 
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7.  Conclusions 

We  have  presented  another  view  of  domain  decomposition,  one  which  may  be  illuminat¬ 
ing  to  those  who  find  the  discrete  presentation  of  the  method  difficult  to  understand.  Using 
this  view,  we  have  analyzed  several  preconditioners  associated  with  the  method.  We  con¬ 
clude,  at  least  for  the  model  problem,  that  the  preconditioners  are  essentially  the  same.  Each 
of  them  contains  some  approximation  to  -{-Ai.^ ,  which  is  the  principle  part  of  the  forward 
operator  of  the  interface  equations.  If  the  aspect  ratios  of  the  component  rectangles  are  not 
excessively  small,  then  each  of  the  preconditioners  can  be  expected  to  do  rather  well. 

We  have  also  examined  a  particular  variable  coefficient  problem  and  introduced  a  new 
preconditioner  suitable  for  it.  One  key  fact  which  can  be  inferred  from  our  analysis  is  that 
piecewise  constant  coefficient  problems  are  readily  amenable  to  solution  using  domain  decom¬ 
position  ideas.  Although  we  examined  a  very  simple  problem,  the  basic  ideas  cany  over  to 
more  complicated  problems.  In  particular,  insight  obtained  from  the  model  problem  may  be 
useful  for  those  implementing  the  algorithm  of  Bramble  Pasciak  and  Hubbard  in  [2]  or  [3]. 

Our  numerical  experiments  demonstrated  that  the  model  problem  is  an  accurate  predic¬ 
tor  of  the  behavior  of  the  method  ”L”  shaped  and  approximately  "L*  shaped  domains.  In 
the  case  of  variable  coefficients,  we  also  confirmed  that  conclusions  based  on  the  model  prob¬ 
lem  were  correct. 
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M2 
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Table  3 

°2 

Condition  number  vs.  — . 

«i 

Qfj  ”  .6  -J —  1.0  ^2  *  1*0 

h  —  .01 

Variable  coefficients,  e  j  =  1.0  cs  =  190.0 


